#!/usr/bin/perl

use Math::Trig;

# Constants
#===============================================================================
$c = 299792458.0;
$deg_to_rad = pi/180.0;

#===============================================================================
# Vermeulen parameters
#$a5   = 0.00382;
#$a6   = 0.00655;
#$p5   = 5.838;
#$p6   = 6.290;
#$t5p8 = 3588.5;
#$t6p3 = 3587.9;

# KMB old best fit parameters 
#$a5   = 0.00439;
#$a6   = 0.00732;
#$t5p8 = 3587.39;
#$t6p3 = 3587.97;
#$p5   = 5.845;
#$p6   = 6.2878; 

# KMB new best fit parameters 
$a5   = 0.00119145971;
$a6   = 0.00517231109;
$t5p8 = 3588.14956;
$t6p3 = 3587.61392;
$p5   = 5.83990251;
$p6   = 6.29175947;


# Advanced inputs
#===============================================================================
# CHOOSE CONVENTION
#   - this is either A or B or C (see attached image)

# CHOOSE IF NODDING IS SWITCHED ON
#   - this nodding is comparable with the stochastic wobbling so sometimes 
#     it's nice to have it turned off

# CHOOSE IF MODEL PARAMETERS DIFFER FROM USUAL
$SPEED = 0.259;
$LOS   = 78.5697;
$PSI   = 19.8420;
$CHI   = 10.0;
$i     = $LOS * $deg_to_rad;
$psi   = $PSI * $deg_to_rad;
$chi   = $CHI * $deg_to_rad; # angle on sky

# [INPUT] Choose date. (Today or date) (data goes back as far as 1978!)
#===============================================================================
$t     = 4500.0;


# Other parameters
#===============================================================================
$s_rot = 1;
$s_jet = 1;
# $s_jet < 0 for right, west jet and $s_jet > 0 for left, east jet 

$nod   = 0;
$v     = $SPEED;
$psi   = $PSI;
$omega = 1.0 * pi; ## ???! WHAT IS THIS
$t_0   = 0;
$tref  = 0;


# CALCULATE & OUTPUT redshift etc
#===============================================================================
$psi5  = 2.0 * pi * ($t - $t5p8)/$p5;
$psi6  = 2.0 * pi * ($t - $t6p3)/$p6;

$beta  = $v;
$gamma = 1.0 / sqrt(1.0 - $beta * $beta);

$p1 = $beta * sin($theta) * sin($i) * cos($phi);
$p2 = $beta * cos($theta) * cos($i);

# NOTE: Where does this equation come from?
$vdotx = $s_jet * $beta * ((sin($psi) * sin($i) * cos($omega*($t_0 - $tref)) + cos($psi) * cos($i))) + $nod * $s_jet * $c * ($a5 * cos($psi5) + $a6 * cos($psi6));

# NOTE: is this equation correct?? where does it come from...?
#$redshift = $gamma - 1.0 - $gamma * $vdotx/$c;
# ... perhaps this is right...?
$redshift = $gamma * (1.0 + $vdotx) - 1.0;

# Printing
#===============================================================================
print "\n";
print "==== Inputs ====\n";
print "* t ...... = $t (MJD?)\n";
print "* t_0 .... = $t_0\n";
print "* tref ... = $tref\n";
print "* v/c .... = $v\n";
print "* c ...... = $c\n";
print "* psi..... = $psi\n";
print "* i ...... = $i\n";
print "* omega .. = $omega\n";

print "\n";
print "==== Results ====\n";
print "* psi5 ... = $psi5\n";
print "* psi6 ... = $psi6\n";
print "* vdotx .. = $vdotx\n";
print "* beta ... = $beta\n";
print "* gamma .. = $gamma\n";
print "* Redshift = $redshift\n";
print "\n";


